% Read raw data file
data = dlmread('jtpa.raw');

% Subset to male observations
male = data(:,5) == 1;
data = data(male,:);

% Define main variables we'll use
y = data(:,2); % earnings
d = dummy(data(:,3)); % randomly assigned offer of job training indicator
maineffect{1} = dummy(data(:,6)); % hs grad
maineffect{2} = dummy(data(:,7)+2*data(:,8));  % 7 = black, 8 = hispanic
maineffect{3} = dummy(data(:,12:16)*((1:5)'));  % age2225, age2629, age3035, age3644, age4554
maineffect{4} = dummy(data(:,9));  % married
maineffect{5} = dummy(data(:,10));  % work less than 13 weeks
maineffect{6} = dummy(data(:,19));  % f2sms??

% Create all the interactions
px = [size(maineffect{1},2) ; size(maineffect{2},2) ; ...
    size(maineffect{3},2) ; size(maineffect{4},2) ; ...
    size(maineffect{5},2) ; size(maineffect{6},2)];

% two-way interactions
cellnum = 0;
for jj1 = 1:6
    for jj2 = (jj1+1):6
        cellnum = cellnum+1;
        tmp = [];
        v1 = maineffect{jj1};
        v2 = maineffect{jj2};
        for kk1 = 1:px(jj1)
            for kk2 = 1:px(jj2)
                tmp = [tmp , v1(:,kk1).*v2(:,kk2)];
            end
        end
        interact2{cellnum} = tmp;
    end
end

% three-way interactions
cellnum = 0;
for jj1 = 1:6
    for jj2 = (jj1+1):6
        for jj3 = (jj2+1):6
            cellnum = cellnum+1;
            tmp = [];
            v1 = maineffect{jj1};
            v2 = maineffect{jj2};
            v3 = maineffect{jj3};
            for kk1 = 1:px(jj1)
                for kk2 = 1:px(jj2)
                    for kk3 = 1:px(jj3) 
                        tmp = [tmp , v1(:,kk1).*v2(:,kk2).*v3(:,kk3)];
                    end
                end
            end
            interact3{cellnum} = tmp;
        end
    end
end

% four-way interactions
cellnum = 0;
for jj1 = 1:6
    for jj2 = (jj1+1):6
        for jj3 = (jj2+1):6
            for jj4 = (jj3+1):6
                cellnum = cellnum+1;
                tmp = [];
                v1 = maineffect{jj1};
                v2 = maineffect{jj2};
                v3 = maineffect{jj3};
                v4 = maineffect{jj4};
                for kk1 = 1:px(jj1)
                    for kk2 = 1:px(jj2)
                        for kk3 = 1:px(jj3)
                            for kk4 = 1:px(jj4)
                                tmp = [tmp , v1(:,kk1).*v2(:,kk2).*v3(:,kk3).*v4(:,kk4)];
                            end
                        end
                    end
                end
                interact4{cellnum} = tmp;
            end
        end
    end
end

% five-way interactions
cellnum = 0;
for jj1 = 1:6
    for jj2 = (jj1+1):6
        for jj3 = (jj2+1):6
            for jj4 = (jj3+1):6
                for jj5 = (jj4+1):6
                    cellnum = cellnum+1;
                    tmp = [];
                    v1 = maineffect{jj1};
                    v2 = maineffect{jj2};
                    v3 = maineffect{jj3};
                    v4 = maineffect{jj4};
                    v5 = maineffect{jj5};
                    for kk1 = 1:px(jj1)
                        for kk2 = 1:px(jj2)
                            for kk3 = 1:px(jj3)
                                for kk4 = 1:px(jj4)
                                    for kk5 = 1:px(jj5)
                                        tmp = [tmp , v1(:,kk1).*v2(:,kk2).*v3(:,kk3).*v4(:,kk4).*v5(:,kk5)];
                                    end
                                end
                            end
                        end
                    end
                    interact5{cellnum} = tmp;
                end
            end
        end
    end
end

% six-way interactions
cellnum = 0;
for jj1 = 1:6
    for jj2 = (jj1+1):6
        for jj3 = (jj2+1):6
            for jj4 = (jj3+1):6
                for jj5 = (jj4+1):6
                    for jj6 = (jj5+1):6
                        cellnum = cellnum+1;
                        tmp = [];
                        v1 = maineffect{jj1};
                        v2 = maineffect{jj2};
                        v3 = maineffect{jj3};
                        v4 = maineffect{jj4};
                        v5 = maineffect{jj5};
                        v6 = maineffect{jj6};
                        for kk1 = 1:px(jj1)
                            for kk2 = 1:px(jj2)
                                for kk3 = 1:px(jj3)
                                    for kk4 = 1:px(jj4)
                                        for kk5 = 1:px(jj5)
                                            for kk6 = 1:px(jj6)
                                                tmp = [tmp , v1(:,kk1).*v2(:,kk2).*v3(:,kk3).*v4(:,kk4).*v5(:,kk5).*v6(:,kk6)];
                                            end
                                        end
                                    end
                                end
                            end
                        end
                        interact6{cellnum} = tmp;
                    end
                end
            end
        end
    end
end

X = [cell2mat(maineffect) , cell2mat(interact2) , cell2mat(interact3) , ...
    cell2mat(interact4) , cell2mat(interact5) , cell2mat(interact6)];

% Drop columns with very observations in either treatment or control state
nx = sum(X(d == 1,:) > 0);
X = X(:,nx > 5);

ndx = sum(X(d == 0,:) > 0);
X = X(:,ndx > 5);

% Drop any columns that would lead to singularity in regressions within
% treatment or control subsets
[~,Rx,ex] = qr(X(d == 1,:),0);
keepvar = ex(abs(diag(Rx)) > 1e-6);
X = X(:,keepvar);

[~,Rx,ex] = qr(X(d == 0,:),0);
keepvar = ex(abs(diag(Rx)) > 1e-6);
X = X(:,keepvar);

% Keep the variables we care about
clearvars -except y d X ;
